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For molecules weakly coupled to leads the exact Kohn-Sham (KS) conductance can be orders of 
magnitude larger than the true conductance due to the lack of dynamical exchange-correlation (xc) 
corrections. In this work we show how to incorporate dynamical effects in KS transport calcula- 
tions. The only quantity needed is the static xc potential in the molecular junction. Our scheme 
provides a comprehensive description of Coulomb blockade without breaking the spin symmetry. 
This is explicitly demonstrated in single- wall nanotubes where the corrected conductance is in good 
agreement with experimental data whereas the KS conductance fails dramatically. 
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The active field of molecular electronics [H remains a 
challenge for ab initio methods. Density Functional The- 
ory (DFT) is at present the only viable route for an atom- 
istic description of complex molecular junctions. How- 
ever, DFT transport calculations still suffer from some 
conceptual difficulties. Especially thorny is the problem 
of weakly coupled molecules where level alignment and 
charging effects play a prominent role. Toher et al. Q 
showed that the discontinuity of the exchange-correlation 
(xc) potential Q is crucial for opening a gap in the I- 
V characteristics of closed-shell molecules, and hence for 
capturing the Coulomb blockade (CB) effect at even elec- 
tron numbers TV. At the same time, the current "DFT 
understanding" of CB in opcn-shcU molecules is rather 
limited: while the discontinuity suppresses the DFT con- 
ductance at even N, it has the opposite effect at odd 
A'^, see below. In this work we provide a comprehensive 
picture of CB, valid for all N without breaking the spin 
symmetry. The key ingredient is the dynamical xc cor- 
rection to the conductance 0, [Ij which, remarkably, can 
be expressed exclusively in terms of static DFT quanti- 
ties. Comparison with experiments on single-wall nan- 
otubes shows a significant improvement of the corrected 
conductances over the DFT results. 

At zero temperature and for single- channel junctions 
the equality Gg = G between the Kohn-Sham (KS) con- 
ductance Gs and the true conductance G is a consequence 
of the Friedel sum rule Q . The key ingredient to repro- 
duce the Kondo plateau in Gs is the discontinuity [7|-l9|. 
However, at temperatures higher than the Kondo tem- 
perature Tk, Gs can be orders of magnitude larger than 
G We can understand this discrepancy by modeling 
the molecule as a single level (HOMO or LUMO) of en- 
ergy V and Coulomb repulsion U coupled to left {L) and 
right [R) featureless leads contributing 7 = Jl+Jb. to the 




FIG. 1. Top: Electron number A'^ versus gate in MB and 
DFT (indistinguishable) for a single level coupled to feature- 
less leads with U = 10, /i = at various temperatures T (all 
energies in units of 7). Bottom: G from Eq. ([2]) (solid) and 
Gb from Eq. Q (dashed) in units of Go = 2e^//i. 



broadening of the spectral peaks. Given the Many-Body 
(MB) spectral function A{u!) the number of electrons is 



N = 2j f{io)A{Lu), 
whereas the zero-bias conductance reads 
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with the Fermi function /(w) = l/(e^('^"^) -I- 1) at in- 
verse temperature /? = 1/T and chemical potential /i. 
For temperatures T ^ Tk the Abrikosov-Suhl (AS) res- 
onance is strongly suppressed and the spectral function 
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is well represented by [lO[ 

A{u) = nL^{uj - u - [/) + (1 - n)L^{u - v) (3) 

where n = N/2 and Lj{u}) = 7/(w^ + 7^/4). Consider 
now the KS system whose spectral function is As(uj) = 
Lj{Ld — V — whxc[^])- The Hartree-xc (Hxc) potential 
I'Hxc is such that the number of electrons N which solves 
N — 2 J f{u})As{uj) is the same as in Eq. ([T|). We obtain 
I'Hxc by reverse engineering. In Fig. [T] (top) we show the 
MB and DFT N-v curves (which arc indistinguishable) 
for three different temperatures. The values of N are 
essentially independent of T up to T ^ 7. Having whxc 
we calculate the KS conductance from 
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In Fig. [T] (bottom) we compare G and Gg ■ Despite the 
fact that the MB and DFT densities are identical there 
is no signature of the CB peaks in Gg- Therefore any 
DFT approach, based either on nonequilibrium Green's 
functions or on scattering formalisms, considerably over- 
estimates the conductance of open-shell molecules. We 
remark that the physical situation discussed here is dis- 
tinct from that of Ref. [H where the discontinuity keeps 
the HOMO doubly occupied and the LUMO empty when 
gating the molecule (closed shell). 

Dynamical xc effects: Open-shell molecules in the CB 
regime are probably the most striking example of the 
failure of static DFT methods. We now derive an ex- 
act formula for G in terms of time-dependent (TD) DFT 
quantities [lli. We take the leads as two jcUia (the ar- 
gument can be generalized to more realistic leads) and 
choose z as the longitudinal coordinate so that z ^ —00 
is in the left lead, a = L, whereas z -> cx) is in the right 
lead, a = R. Let SV^ be the variation in the potential 
of lead a. This perturbation generates a current [1, [Hj 
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where SV^^^ = limt_^oo hm^^j^oo '5fHxc(r, i), sr/l = ±, 
is the asymptotic value of the variation of the Hxc po- 
tential (Swhxc- From linear-response TDDFT 



lim lim 



dr'dt'fH.cir,r';t^t')Sn{r',t') 

(6) 

where /hxc is the TDDFT kernel and 6n{r,t) is the den- 
sity variation. The assumption of a steady state im- 
plies that /hxc for |t — i'| — > 00 and Sn{r,t 
00) = SaSn for r in lead a. In Eq. ([6]) the contri- 
bution of the molecular region to the spatial integral 
is negligible in the thermodynamic limit. If we dc- 
fii^c /hxc = liint->oo lim^^s^oo /i^ad /3 '^i"'/Hxc(r, r'; t) then 



ilcad 13 

'Hxc - E/3=L,i?,/Hxc'S/3^"- We emphasize that /^^^ is 
not the static DFT kernel since the limit i — > 00 is taken 
after the limit |z| — 00 and these two limits, in general. 



do not commute [l2| . This implies that we cannot model 



/hxc performing DFT calculations on leads of finite 
length. Inserting the expression for SV^^^ into Eq. ([5]) 
we find SI = {SV^ - SV")Gs - ^GsSn where 
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The expression for 51 is correctly gauge invariant. The 
kernel /hxc is defined up to the addition of an arbitrary 
function ^(r) + g{v') [l3| and we see that $ is invariant 
under this transformation. In conclusion 



G = 



SI 
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The quantity x = Sn/SI ^ 1/{vf<^) with vp the Fermi 



IJ. Thus 



velocity and a the cross section of the leads 
all dynamical xc effects are contained in $. 

Approximations to $: By definition the product x$ = 
((^^Hxc - '5V[f^J/((5F^ - 5V^) is the ratio between the 
Hxc bias and the physical bias. To gain some insight 
into its density dependence we consider again the single 
level model. For N ^ 1 the real and KS systems respond 
similarly and consequently G ~ G^. On the other hand 
for = 1 we have G ~ whereas Gg ~ Go = 2e^//i 
the quantum of conductance. Therefore $ is small for 
N 1 and large for = 1. Interestingly the quantity 
dvuxc/dN behaves similarly. Is there any relation be- 
tween $ and dvB_xc/dN7 If so this relation would simplify 
enormously the problem of estimating the dynamical xc 
correction since dvB_xc/dN can be calculated from static 
DFT. In the following we show that in the CB regime 
this relation does actually exist. 

Consider the system in equilibrium. Using Eq. 
the compressibility k = dN/dfi can be written as k = 



rG + 2//(o 



,dA(a 



, where we identified the conduc- 
If we define the quantity R = 
G-r^. As the MB and 
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tance G of Eq. ^ 
-2//(.)^then. = ^ 
DFT densities are the same, the MB and DFT com- 



-G. 



pressibilities are the same too. Hence k 

2 J f{u})^^^gj^, where we identified the KS conductance 
Gs of Eq. (U]). The KS spectral function depends on /i 
through N, and the dependence on N is all contained 
in «Hxc: ^^ = Using this resuh under 

the integral sign, solving for k and equating the MB and 
DFT expressions we get 
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No approximations have been made so far. Let us study 
the dependence of R on temperature. 

We first consider the low temperature case. For sim- 
plicity we take 7^ = and set v = —U/2 at the 
particle-hole (ph) symmetric point. At zero tempera- 
ture G = G,. = Go and hence R = ■ The 
correctness of this equality can easily be verified. The 



compressibility is weakly dependent on temperature up 
to r ~ 7 [lij. Therefore the uhxc obtained by reverse 
engineering is a good approximation even for T = 0. 
At the ph symmetric point {N = 1) this approxima- 
tion gives dv-a^c/dN ~ C/V7 0. At T = the 
MB spectral function for lo ^ fi is dominated by the 
AS resonance A{u! • 



/i) ~ —^LT^i^ — Ai) and thus 
For the compressibility we have 
from the Bethe-Ansatz [H k = {8"f)/{'KU^) [1 + 0{j/U)] 
from which it follows that also R goes like (C//7)^. 



For 



temperatures T > Tk the AS resonance broadens and its 
height decreases like /i(T/Tk) where his a universal func- 
tion which approaches zero at high T [T^]- This means 
that R - h{T/TK){UhY remains large until the AS 
resonance disappears. No simple relation between $ and 
dv-nxc/dN exists when Kondo correlations arc present. 

At temperatures T ^ Tk thermal fluctuations destroy 
the Kondo effect and the MB spectral function is well ap- 
proximated by Eq. ©. Therefore R{v) = I{v)-I[v + U) 
where I{E) = ] f{uj)L^{uj ~ E). Hence R{v) - 
for u > ^ (iV - 0) or w < ^ - [/ (A^ - 2) and 
R{v) < 1 otherwise. We can derive a more conve- 
nient expression for R by inserting Eq. ^ into Eq. 
(HD to find N ^ N liv + U) + {2 - N)I{v) and hence 
1 + R = 2I{v)/N. Unfortunately I{v) is not an ex- 
plicit "density functional" due to the implicit depen- 
dence oi V = v[N]. However, for TV < 1 a good ap- 
proximation to is ~ 2/(7j)/(l + I{v)) which implies 
I{v) = N/ (2 — A^) . Taking into account the ph symmetry 
we can approximate R by an explicit density functional 
as 1 + i? 2/(1 + \6N\) where 6N = N - 1. Inserting 
this into Eq. ^ we deduce our main result 
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which provides a simple correction to the KS conduc- 
tance. Remarkably, the dynamical xc correction of 
Eq. ([5]) is entirely expressed in terms of static DFT quan- 
tities of the molecular junction. In Fig. [5] we compare 
the conductance G of Eq. Q with the conductance cal- 
culated from Eq. (fTO)l . Even though the approximate R 
is not on top of the exact one, see inset, the agreement 
between the two conductances is extremely good. The 
position, width and height of the peaks as well as the 
decay for large |?;| are all well reproduced. In fact, the 
approximate R is rather accurate except in the region 
where A^ 1 (/i — C/ < w < /i). In this region, how- 
ever, dvuxc/dN is large and the small error in R docs 
not significantly affect G. 

Application to physical systems: In real molecules 
whxc is a r-dependent functional of the density. We 
write t'Hxc(r) = Svuxd^) + vhxc as the sum of a func- 
tional (5uhxc with a weak dependence on N ~ Jy drn(r) 
and a spatially uniform part whxc — y Jy dr vn^ci^) , 
where the integral is over the volume V of the molec- 
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FIG. 2. Conductance from Eq. Q (solid) and from Eq. (fTU)) 
(dashed). The inset shows a comparison between the exact 
and the approximate R. Same parameters as in Fig. [T] 



ular junction. For weakly coupled molecules whxc ex- 
hibits sharp steps as a function A^ when N crosses an 
integer. If we are at resonance and spin fluctuations 
are suppressed (no Kondo effect) then the KS conduc- 
tance should be corrected according to Eq. pO[) in which 
dv-nxc/dN —5- dvxc/dN. Below wc argue that this cor- 
rection applies out-of-rcsonance too. Let fi be in the 
HOMO-LUMO gap and consider a two-level system with 
Fq the 2x2 broadening matrix. For general Fq, no sim- 
ple analytic relation between G and A^ exists. However 
if r„,,„, = {-fa/2)5mi then iV = 2 / /(w)Tr [A(w)] and 
G = -2^^ / /'(w)Tr[A(w)]. Discarding the depen- 
dence of (5vhxc on fi (which is weak by definition) wc can 
go through the same steps of the single-level derivation 
and find again Eq. (fTU|l . It is therefore reasonable to ex- 
pect that the KS conductance should be corrected even 
out-of-resonance (closed shell) and that this correction 
should be proportional to Ggdvuxc/dN . 

We propose a scheme to calculate G from DFT. Given 
the KS Hamiltonian matrix /iKS.mZ = SmiQ and the 
broadening matrices Ta.mi we determine the density and 
Gs in the usual manner. G is then obtained from Eq. (jlOp 
where SN is the deviation of (A — Int[A^]) from 1 whereas 
7a = 7a(-A) ~ Ta,HH H jJL ^ eR (rcsonancc, open shell) 
and 7a (A) ~ \{Ta,HH + ^a^LL) if /.t ~ ^(cL + en) is in 
the HOMO-LUMO gap (out-of-resonance, closed sheU). 
One could improve the approximation to 7a using differ- 
ent weights, but the qualitative features of the results are 
independent of these details, see below. 

To illustrate the importance of the dynamical xc cor- 
rection we consider two paradigmatic junctions in which 
5v}ixc can be discarded. For ?;hxc we choose a best fit of 
the zero-temperature limit of the single-level Hxc poten- 
tial, but now sum over all possible charged states of the 




FIG. 3. KS conductance (dashed) and TDDFT conductance 
from Eq. ^ (solid) for the HOMO-LUMO model with diag- 
onal and off-diagonal F-matrices (left axis) and KS energies 
^H/L = eoH/L + VHxc (right axis). 



molecule [l9[, i.e., 



v-^C/(X) fN-K\ 

> arctan ; — - 

^ IT \W{K)) 



(11) 



The charging energies U{N) are given by the xc part 
of the derivative discontinuity of the molecule with 
N electrons [l]. For the widths we take W{N) = 
0.167(A'^)/;7(7V) which is consistent with Ref.[ll 

HOMO-L UMO model: We consider a two- level system 
which has 2 electrons in the HOMO in the charge neu- 
tral state. Let eon ~ ~eoL = ~eo < be the nonin- 

teracting single-particle energies, = 2 \ \ 

and U{N) = U independent of N. We solve the self- 
consistent equation for the density with U — 10, eq — 5, 
/i = and (3 — 10 (all energies in units of 7). Gg and G 
from Eq. (fTO|) are shown in Fig. [3] (left axis). As expected 
the discontinuity of vuxc opens a gap in Gs for even N, 
in agreement with the results of Ref. 3- Here the dynam- 
ical xc correction only weakly affects Gs since dvny^c/dN 
is multiplied by "C 1. On the other hand for odd N 
the KS conductance exhibits a Kondo plateau due to the 
pinning of the KS level to /i, see right axis. This is the 
regime previously discussed and no CB is observed. The 
dynamical xc correction remedies this deficiency. We re- 
mark that the results remain essentially unaltered if the 
off-diagonal matrix elements of Ta are discarded which 
amounts to neglecting interference effects pol |. 

SWNT: Experimental evidence of CB oscillations has 
been recently reported in metallic single- wall nanotubes 
(SWNT) quantum dots [2l|,[2l|. The finite length of the 
SWNT causes a level quantization of the twofold degen- 
erate bands. Following Ref. |2l| we model the noninter- 
acting single-particle energy levels as eoi,y = lA + l{i'~l)S 
where 1/ = 0, 1 is the band index and / is the integer of 




FIG. 4. KS and TDDFT conductance (Eq. dTO])) for a SWNT 
quantum dot in comparison to experimental conductance 
from Ref. [2ll . as function of gate voltage Vg . 



the quantized wavevector along the nanotube axis. Since 
the wavevector is a good quantum number our approxi- 
mation (5whxc = is justified. The constant interaction 
model 



241 to ac- 



[23| has been refined by Oreg et al. 
count for the observed fourfold periodicity in the elec- 
tron addition energy. We choose the U{Kys in Eq. pT|) 
according to 0: U[l) = U{3) = Ec + SU + J, and 
U{2) = U{A) = Ec- SU, and, for K > 4, U{K) = U{i) 
li K — i ~0 mod4 with z = 1, 2, 3, 4. The average values 
of these parameters can be found in Ref. [2TI. Here we 
change them slightly to match the peak positions of the 
SWNT of length ~ 100 nm (all energies are in meV): 
A = 9.2, 5 = 2.27, charging energy Ec = 2.485, ex- 
change energy J = 0.7, extra charging energy for doubly 
occupied levels 5U = 0.37. For the broadening matrix 
we make the approximation Ta^mi ~ {l /'2)Smi (no visible 
interference from the experiment). In Fig. 2] we compare 
the KS, TDDFT and experimental conductance versus 
the gate voltage Vg = vq + av where vq is chosen in or- 
der to have the same reference energy as in Ref. |2l| and 
we estimated the ratio a — C/Gg ^ 250 between the 
total capacitance and the gate capacitance from the ex- 
perimental data. In our calculation we have taken into 
account 22 single-particle levels and chose a temperature 
T ^ 7 ^ 1 meV. Despite the simplicity of the model 
we clearly see that the dynamical xc correction improves 
considerably over Gs which exhibits two deformed Kondo 
plateaus per period and misses the fourfould periodicity. 
The qualitative behavior of G and Gs does not change by 
varying the parameters within a reasonable range around 
the average values reported in Ref. [2l| (not shown). 

In conclusion we highlighted the importance of the dis- 
continuity for an accurate Gs and for an accurate dynam- 
ical xc correction to G. Approaches to generate discon- 
tinuous xc potentials are emerging both in the static [2^ 
and dynamical [23| case. Our scheme provides a coher- 
ent picture of CB within (TD)DFT without breaking the 
spin symmetry. By application to two different model 
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molecular junctions we showed that the dynamical xc 
correction reduces Gg thus contributing to close the gap 
between theory and experiments. 
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